Load data

pcv7_st <- c("4","6B","9V","14","18C","19F","23F")
pcv13_st <- c("1","3","4","5","6A","6B","7F","9V","14","18C","19A","19F","23F")
pcv15_st <- c("1","3","4","5","6A","6B","7F","9V","14","18C","19A","19F","22F","23F","33F")
pcv20_st <- c("1","3","4","5","6A","6B","7F","8","9V","10A","11A","12F","14","15B","18C","19A","19F","22F","23F","33F")

data_raw <- read_rds("data/ABCs_st_1998_2021.rds")

data <- data_raw |> 
  rename(agec = "Age.Group..years.",
         year = Year,
         st = IPD.Serotype,
         N_IPD = Frequency.Count)  |> 
  mutate(st = if_else(st == '16', '16F', st))  |> 
  group_by(st, year) |> 
  summarize(N_IPD = sum(N_IPD))  |> 
  ungroup() |> 
  complete(year, st, fill = list(N_IPD = 0)) |> 
  mutate(pcv7 = if_else(st %in% pcv7_st, 1, 0),
         pcv13 = if_else(st %in% pcv13_st, 1, 0),
         pcv15 = if_else(st %in% pcv15_st, 1, 0),
         pcv20 = if_else(st %in% pcv20_st, 1, 0),
         firstvax = case_when(pcv7 == 1 ~ "pcv7",
                              pcv13 == 1 ~ "pcv13",
                              pcv15 == 1 ~ "pcv15",
                              pcv20 == 1 ~ "pcv20",
                              .default = ""))
## `summarise()` has grouped output by 'st'. You can override using the `.groups`
## argument.

1: Only PCV-7 serotypes

Wrangle data

# wrangle dataframe
mod1_df <- data |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  filter(st %in% pcv7_st)  |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

# convert to matrix
mod1_mat <- mod1_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

Define model string

mod1_string <- "
  model {
    
    for(i in 1:N_years){
      
      for(j in 1:N_sts){
        N_IPD[i,j] ~ dnegbin(prob[i,j], r[j])
        
        prob[i,j]<- r[j] / (r[j] + lambda[i,j])  # Likelihood 
        
        log(lambda[i,j]) <- epsilon1[i,j]        # serotype-specific intercept + AR(1) effect centered around global effects
        }
      }
    
    # Global AR(1) effect
        beta1[1] ~ dnorm(alpha1, (1 - rho_beta^2) * tau_beta1)
        
    for(i in 2:N_years){
      beta1[i] ~ dnorm(alpha1 + rho_beta * beta1[i-1], tau_beta1)
      }
        
    # Serotype-specific AR(1) effect
    
    for(j in 1:N_sts){
      epsilon1[1,j] ~ dnorm(delta1[j] + beta1[1], (1 - rho_eps^2) * tau_eps)
      
      for(i in 2:N_years){
        epsilon1[i,j] ~ dnorm(delta1[j] + beta1[i] + rho_eps * epsilon1[i-1, j], tau_eps)
        }
      }
        
    # Priors
    alpha1 ~ dnorm(0, 1e-4)
    tau_global ~ dgamma(0.01,0.01)

    rho_beta ~ dunif(-1, 1)   # Uniform prior for rho_beta -- global AR(1)
    rho_eps ~ dunif(-1, 1)    # Prior for rho_eps same for all STs
    
    tau_beta1 ~ dgamma(3, 2)  # Tight prior for tau
    tau_eps ~ dgamma(3, 2)    # Tight prior for tau, shared for all serotypes
    
    for(j in 1:N_sts){
       delta1[j] ~ dnorm(0, tau_global)  # serotype means centered around 0
      
       r[j] ~ dunif(0, 250) #serotype dispersion parameter
      }
}"

Specify initial values

inits1 = list(".RNG.seed"=c(123), ".RNG.name"='base::Wichmann-Hill')
inits2 = list(".RNG.seed"=c(456), ".RNG.name"='base::Wichmann-Hill')
inits3 = list(".RNG.seed"=c(789), ".RNG.name"='base::Wichmann-Hill')

Create JAGS model object

mod1 <- jags.model(textConnection(mod1_string),
                   data = list("N_IPD" = mod1_mat,
                               "N_years" = max(mod1_df$year_n),
                               "N_sts" = ncol(mod1_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

params1 <- c('alpha1', 'epsilon1', 'delta1', 'beta1', 'tau_beta1', 'lambda')

mod1_postsamp <- coda.samples(mod1, params1, n.iter = 10000)

# save(mod1_postsamp, file = "data/interim/mod1_postsamp.rda")

Plot obs vs pred

load("data/interim/mod1_postsamp.rda")

st_mapping <- setNames(colnames(mod1_mat), 1:length(pcv7_st))
year_mapping <- setNames(unique(mod1_df$year), 1:25)

mod1_summary <- gather_draws(mod1_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping[as.character(st)],
         year = year_mapping[as.character(year)]) |> 
  left_join(mod1_df, by = c("year", "st"))
mod1_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod1_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod1_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod1_postsamp, ask = TRUE)

Plot trajectory of global AR(1)

mod1_beta1 <- gather_draws(mod1_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping[as.character(i)])

mod1_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

2: PCV-7 serotypes plus five uncommon serotypes

data |> 
  filter(!(st %in% pcv7_st)) |> 
  group_by(st) |> 
  summarise(N_IPD = sum(N_IPD)) |> 
  arrange(N_IPD)
## # A tibble: 80 × 2
##    st    N_IPD
##    <chr> <int>
##  1 12B       1
##  2 15D       1
##  3 17        1
##  4 17A       1
##  5 18        1
##  6 19B       1
##  7 19C       1
##  8 25F       1
##  9 30        1
## 10 33B       1
## # ℹ 70 more rows
data |> 
  complete(year, st, fill = list(N_IPD = 0)) |> 
  ggplot(aes(x = year, y = N_IPD)) +
  geom_line() +
  facet_wrap(~st) +
  theme_bw()

Wrangle data

# define  serotypes to include in model
pcv7plus_st <- c("4","6B","9V","14","18C","19F","23F","15F","33A","35A","9L","29")

# wrangle dataframe
mod2_df <- data |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  filter(st %in% pcv7plus_st)  |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

# convert to matrix
mod2_mat <- mod2_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

Create JAGS model object

mod2 <- jags.model(textConnection(mod1_string),
                   data = list("N_IPD" = mod2_mat,
                               "N_years" = max(mod2_df$year_n),
                               "N_sts" = ncol(mod2_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod2_postsamp <- coda.samples(mod2, params1, n.iter = 10000)

# save(mod2_postsamp, file = "data/interim/mod2_postsamp.rda")

Plot obs vs pred

load("data/interim/mod2_postsamp.rda")

st_mapping2 <- setNames(colnames(mod2_mat), 1:length(pcv7plus_st))
year_mapping2 <- setNames(unique(mod2_df$year), 1:25)

mod2_summary <- gather_draws(mod2_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping2[as.character(st)],
         year = year_mapping2[as.character(year)]) |> 
  left_join(mod2_df, by = c("year", "st"))
mod2_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod2_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod2_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod2_postsamp, ask = TRUE)
# traceplot(mod2_postsamp[, grep("^lambda", varnames(mod2_postsamp))])

Plot trajectory of global AR(1)

mod2_beta1 <- gather_draws(mod2_postsamp, beta1[i]) |>
  median_hdi() |>
  mutate(year = year_mapping2[as.character(i)])

mod2_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

3: Most-common non-PCV-7 serotypes

Wrangle data

# define  serotypes to include in model
nonpcv7_st <- c("3","19A","22F","7F","12F",
                "9N","33F","11A","23A","16F",
                "6C","35B","8","15A","6A",
                "23B","20","31","15B","10A")

# wrangle dataframe
mod3_df <- data |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  filter(st %in% nonpcv7_st)  |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

# convert to matrix
mod3_mat <- mod3_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

Create JAGS model object

mod3 <- jags.model(textConnection(mod1_string),
                   data = list("N_IPD" = mod3_mat,
                               "N_years" = max(mod3_df$year_n),
                               "N_sts" = ncol(mod3_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod3_postsamp <- coda.samples(mod3, params1, n.iter = 10000)

# save(mod3_postsamp, file = "data/interim/mod3_postsamp.rda")

Plot obs vs pred

load("data/interim/mod3_postsamp.rda")

st_mapping3 <- setNames(colnames(mod3_mat), 1:length(nonpcv7_st))
year_mapping3 <- setNames(unique(mod3_df$year), 1:25)

mod3_summary <- gather_draws(mod3_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping3[as.character(st)],
         year = year_mapping3[as.character(year)]) |> 
  left_join(mod3_df, by = c("year", "st"))
mod3_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod3_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod3_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod3_postsamp, ask = TRUE)
# traceplot(mod3_postsamp[, grep("^lambda", varnames(mod3_postsamp))])

Plot trajectory of global AR(1)

mod3_beta1 <- gather_draws(mod3_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping3[as.character(i)])

mod3_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

4: Most-common non-PCV-7 serotypes plus five uncommon serotypes

Wrangle data

# define  serotypes to include in model
nonpcv7plus_st <- c("3","19A","22F","7F","12F",
                    "9N","33F","11A","23A","16F",
                    "6C","35B","8","15A","6A",
                    "23B","20","31","15B","10A",
                    "15F","33A","35A","9L","29")

# wrangle dataframe
mod4_df <- data |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  filter(st %in% nonpcv7plus_st)  |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

# convert to matrix
mod4_mat <- mod4_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

Create JAGS model object

mod4 <- jags.model(textConnection(mod1_string),
                   data = list("N_IPD" = mod4_mat,
                               "N_years" = max(mod4_df$year_n),
                               "N_sts" = ncol(mod4_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod4_postsamp <- coda.samples(mod4, params1, n.iter = 10000)

# save(mod4_postsamp, file = "data/interim/mod4_postsamp.rda")

Plot obs vs pred

load("data/interim/mod4_postsamp.rda")

st_mapping4 <- setNames(colnames(mod4_mat), 1:length(nonpcv7plus_st))
year_mapping4 <- setNames(unique(mod4_df$year), 1:25)

mod4_summary <- gather_draws(mod4_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping4[as.character(st)],
         year = year_mapping4[as.character(year)]) |> 
  left_join(mod4_df, by = c("year", "st"))
mod4_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod4_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod4_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod4_postsamp, ask = TRUE)
# traceplot(mod4_postsamp[, grep("^lambda", varnames(mod4_postsamp))])

Plot trajectory of global AR(1)

mod4_beta1 <- gather_draws(mod4_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping4[as.character(i)])

mod4_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

5: Changing parameters of AR1 - increased precision for tau_eps (Most-common non-PCV-7 serotypes plus five uncommon serotypes)

mod2_string <- "
  model {
    
    for(i in 1:N_years){
      
      for(j in 1:N_sts){
        N_IPD[i,j] ~ dnegbin(prob[i,j], r[j])
        
        prob[i,j]<- r[j] / (r[j] + lambda[i,j])  # Likelihood 
        
        log(lambda[i,j]) <- epsilon1[i,j]        # serotype-specific intercept + AR(1) effect centered around global effects
        }
      }
    
    # Global AR(1) effect
        beta1[1] ~ dnorm(alpha1, (1 - rho_beta^2) * tau_beta1)
        
    for(i in 2:N_years){
      beta1[i] ~ dnorm(alpha1 + rho_beta * beta1[i-1], tau_beta1)
      }
        
    # Serotype-specific AR(1) effect
    
    for(j in 1:N_sts){
      epsilon1[1,j] ~ dnorm(delta1[j] + beta1[1], (1 - rho_eps^2) * tau_eps)
      
      for(i in 2:N_years){
        epsilon1[i,j] ~ dnorm(delta1[j] + beta1[i] + rho_eps * epsilon1[i-1, j], tau_eps)
        }
      }
        
    # Priors
    alpha1 ~ dnorm(0, 1e-4)
    tau_global ~ dgamma(0.01,0.01)

    rho_beta ~ dunif(-1, 1)   # Uniform prior for rho_beta -- global AR(1)
    rho_eps ~ dunif(-1, 1)    # Prior for rho_eps same for all STs
    
    tau_beta1 ~ dgamma(3, 2)  # Tight prior for tau
    tau_eps ~ dgamma(10, 5)   # INCREASED PRECISION FOR TAU, shared for all serotypes
    
    for(j in 1:N_sts){
       delta1[j] ~ dnorm(0, tau_global)  # serotype means centered around 0
      
       r[j] ~ dunif(0, 250) #serotype dispersion parameter
      }
}"

Create JAGS model object

mod5 <- jags.model(textConnection(mod2_string),
                   data = list("N_IPD" = mod4_mat,
                               "N_years" = max(mod4_df$year_n),
                               "N_sts" = ncol(mod4_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod5_postsamp <- coda.samples(mod5, params1, n.iter = 10000)

# save(mod5_postsamp, file = "data/interim/mod5_postsamp.rda")

Plot obs vs pred

load("data/interim/mod5_postsamp.rda")

mod5_summary <- gather_draws(mod5_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping4[as.character(st)],
         year = year_mapping4[as.character(year)]) |> 
  left_join(mod4_df, by = c("year", "st"))
mod5_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod5_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod5_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod4_postsamp, ask = TRUE)
# traceplot(mod4_postsamp[, grep("^lambda", varnames(mod4_postsamp))])

Plot trajectory of global AR(1)

mod5_beta1 <- gather_draws(mod5_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping4[as.character(i)])

mod5_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

6: Changing parameters of AR1 - decreased precision for tau_eps (Most-common non-PCV-7 serotypes plus five uncommon serotypes)

mod3_string <- "
  model {
    
    for(i in 1:N_years){
      
      for(j in 1:N_sts){
        N_IPD[i,j] ~ dnegbin(prob[i,j], r[j])
        
        prob[i,j]<- r[j] / (r[j] + lambda[i,j])  # Likelihood 
        
        log(lambda[i,j]) <- epsilon1[i,j]        # serotype-specific intercept + AR(1) effect centered around global effects
        }
      }
    
    # Global AR(1) effect
        beta1[1] ~ dnorm(alpha1, (1 - rho_beta^2) * tau_beta1)
        
    for(i in 2:N_years){
      beta1[i] ~ dnorm(alpha1 + rho_beta * beta1[i-1], tau_beta1)
      }
        
    # Serotype-specific AR(1) effect
    
    for(j in 1:N_sts){
      epsilon1[1,j] ~ dnorm(delta1[j] + beta1[1], (1 - rho_eps^2) * tau_eps)
      
      for(i in 2:N_years){
        epsilon1[i,j] ~ dnorm(delta1[j] + beta1[i] + rho_eps * epsilon1[i-1, j], tau_eps)
        }
      }
        
    # Priors
    alpha1 ~ dnorm(0, 1e-4)
    tau_global ~ dgamma(0.01,0.01)

    rho_beta ~ dunif(-1, 1)   # Uniform prior for rho_beta -- global AR(1)
    rho_eps ~ dunif(-1, 1)    # Prior for rho_eps same for all STs
    
    tau_beta1 ~ dgamma(3, 2)  # Tight prior for tau
    tau_eps ~ dgamma(1, 0.1)    # DECREASED PRECISION FOR TAU, shared for all serotypes
    
    for(j in 1:N_sts){
       delta1[j] ~ dnorm(0, tau_global)  # serotype means centered around 0
      
       r[j] ~ dunif(0, 250) #serotype dispersion parameter
      }
}"

Create JAGS model object

mod6 <- jags.model(textConnection(mod3_string),
                   data = list("N_IPD" = mod4_mat,
                               "N_years" = max(mod4_df$year_n),
                               "N_sts" = ncol(mod4_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod6_postsamp <- coda.samples(mod6, params1, n.iter = 10000)

# save(mod6_postsamp, file = "data/interim/mod6_postsamp.rda")

Plot obs vs pred

load("data/interim/mod6_postsamp.rda")

mod6_summary <- gather_draws(mod6_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping4[as.character(st)],
         year = year_mapping4[as.character(year)]) |> 
  left_join(mod4_df, by = c("year", "st"))
mod6_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod6_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod6_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

Plot trajectory of global AR(1)

mod6_beta1 <- gather_draws(mod6_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping4[as.character(i)])

mod6_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

7: Changing parameters of AR1 - reducing range of rho_eps (Most-common non-PCV-7 serotypes plus five uncommon serotypes)

mod4_string <- "
  model {
    
    for(i in 1:N_years){
      
      for(j in 1:N_sts){
        N_IPD[i,j] ~ dnegbin(prob[i,j], r[j])
        
        prob[i,j]<- r[j] / (r[j] + lambda[i,j])  # Likelihood 
        
        log(lambda[i,j]) <- epsilon1[i,j]        # serotype-specific intercept + AR(1) effect centered around global effects
        }
      }
    
    # Global AR(1) effect
        beta1[1] ~ dnorm(alpha1, (1 - rho_beta^2) * tau_beta1)
        
    for(i in 2:N_years){
      beta1[i] ~ dnorm(alpha1 + rho_beta * beta1[i-1], tau_beta1)
      }
        
    # Serotype-specific AR(1) effect
    
    for(j in 1:N_sts){
      epsilon1[1,j] ~ dnorm(delta1[j] + beta1[1], (1 - rho_eps^2) * tau_eps)
      
      for(i in 2:N_years){
        epsilon1[i,j] ~ dnorm(delta1[j] + beta1[i] + rho_eps * epsilon1[i-1, j], tau_eps)
        }
      }
        
    # Priors
    alpha1 ~ dnorm(0, 1e-4)
    tau_global ~ dgamma(0.01,0.01)

    rho_beta ~ dunif(-1, 1)   # Uniform prior for rho_beta -- global AR(1)
    rho_eps ~ dunif(-0.5, 0.5)    # Prior for rho_eps same for all STs
    
    tau_beta1 ~ dgamma(3, 2)  # Tight prior for tau
    tau_eps ~ dgamma(3, 2)    # Tight prior for tau, shared for all serotypes
    
    for(j in 1:N_sts){
       delta1[j] ~ dnorm(0, tau_global)  # serotype means centered around 0
      
       r[j] ~ dunif(0, 250) #serotype dispersion parameter
      }
}"

Create JAGS model object

mod7 <- jags.model(textConnection(mod4_string),
                   data = list("N_IPD" = mod4_mat,
                               "N_years" = max(mod4_df$year_n),
                               "N_sts" = ncol(mod4_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod7_postsamp <- coda.samples(mod7, params1, n.iter = 10000)

# save(mod7_postsamp, file = "data/interim/mod7_postsamp.rda")

Plot obs vs pred

load("data/interim/mod7_postsamp.rda")

mod7_summary <- gather_draws(mod7_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping4[as.character(st)],
         year = year_mapping4[as.character(year)]) |> 
  left_join(mod4_df, by = c("year", "st"))
mod7_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(firstvax ~ st, ncol = 5, axes = "all_x") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

mod7_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(firstvax ~ st, ncol = 5, axes = "all_x") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

mod7_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(firstvax ~ st, scales = "free_y", ncol = 5, axes = "all_x") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

Plot trajectory of global AR(1)

mod7_beta1 <- gather_draws(mod7_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping4[as.character(i)])

mod7_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

# run in console
# plot(mod4_postsamp, ask = TRUE)
# traceplot(mod4_postsamp[, grep("^lambda", varnames(mod4_postsamp))])

8: Changing parameters of AR1 - reducing range of rho_eps (PCV-7 serotypes plus five uncommon serotypes)

(did not converge)

Create JAGS model object

mod8 <- jags.model(textConnection(mod4_string),
                   data = list("N_IPD" = mod2_mat,
                               "N_years" = max(mod2_df$year_n),
                               "N_sts" = ncol(mod2_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod8_postsamp <- coda.samples(mod8, params1, n.iter = 50000)

# save(mod8_postsamp, file = "data/interim/mod8_postsamp.rda")
load("data/interim/mod8_postsamp.rda")

beta_vars8 <- grep("^beta", varnames(mod8_postsamp), value = TRUE)
for (i in seq_along(beta_vars8)) {
  cat(paste0("#### ",i," \n"))
  traceplot(mod8_postsamp[, beta_vars8[i]])
  cat("```\n")
}

1

#### 2 <img src="analysis2_files/figure-html/unnamed-chunk-62-2.png" width="672" /> #### 3 #### 4 <img src="analysis2_files/figure-html/unnamed-chunk-62-4.png" width="672" /> #### 5 #### 6 <img src="analysis2_files/figure-html/unnamed-chunk-62-6.png" width="672" /> #### 7 #### 8 <img src="analysis2_files/figure-html/unnamed-chunk-62-8.png" width="672" /> #### 9 #### 10 <img src="analysis2_files/figure-html/unnamed-chunk-62-10.png" width="672" /> #### 11 #### 12 <img src="analysis2_files/figure-html/unnamed-chunk-62-12.png" width="672" /> #### 13 #### 14 <img src="analysis2_files/figure-html/unnamed-chunk-62-14.png" width="672" /> #### 15 #### 16 <img src="analysis2_files/figure-html/unnamed-chunk-62-16.png" width="672" /> #### 17 #### 18 <img src="analysis2_files/figure-html/unnamed-chunk-62-18.png" width="672" /> #### 19 #### 20 <img src="analysis2_files/figure-html/unnamed-chunk-62-20.png" width="672" /> #### 21 #### 22 <img src="analysis2_files/figure-html/unnamed-chunk-62-22.png" width="672" /> #### 23 #### 24 <img src="analysis2_files/figure-html/unnamed-chunk-62-24.png" width="672" /> #### 25 ```

Plot obs vs pred

mod8_summary <- gather_draws(mod8_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping2[as.character(st)],
         year = year_mapping2[as.character(year)]) |> 
  left_join(mod2_df, by = c("year", "st"))

A

mod8_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

B

mod8_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

C

mod8_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod4_postsamp, ask = TRUE)
# traceplot(mod4_postsamp[, grep("^lambda", varnames(mod4_postsamp))])

Plot trajectory of global AR(1)

mod8_beta1 <- gather_draws(mod8_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping2[as.character(i)])

mod8_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1)") +
  theme_minimal()

9: Mixture model with PCV-13

(did not converge)

Wrangle data

# wrangle dataframe
mod9_df <- data |> 
  complete(st, year, fill = list(N_IPD = 0))  |> 
  arrange(st, year) |> 
  filter(st %in% pcv13_st)  |> 
  group_by(st) |> 
  arrange(year) |> 
  mutate(year_n = row_number()) |> 
  ungroup()

# convert to matrix
mod9_mat <- mod9_df |> 
  reshape2::dcast(year ~ st, value.var = 'N_IPD') |> 
  dplyr::select(-year)

Define model string

mod9_string <- "
  model {
    
    for(i in 1:N_years){
      
      for(j in 1:N_sts){
        N_IPD[i,j] ~ dnegbin(prob[i,j], r[j])
        
        prob[i,j] <- r[j] / (r[j] + lambda[i,j])  # Likelihood 
        
        log(lambda[i,j]) <- epsilon1[i,j]  # Serotype-specific intercept + AR(1) effect
      }
    }
    
    # Global AR(1) effect for beta1
    
    beta1[1] ~ dnorm(alpha1, (1 - rho_beta1^2) * tau_beta1)
    for(i in 2:N_years){
      beta1[i] ~ dnorm(alpha1 + rho_beta1 * beta1[i-1], tau_beta1)
    }
    
    # Global AR(1) effect for beta2
    
    beta2[1] ~ dnorm(alpha2, (1 - rho_beta2^2) * tau_beta2)
    for(i in 2:N_years){
      beta2[i] ~ dnorm(alpha2 + rho_beta2 * beta2[i-1], tau_beta2)
    }
        
    # Serotype-specific AR(1) effect with group selection
    for(j in 1:N_sts){
      epsilon1[1,j] ~ dnorm(delta1[j] + beta1[1] * grp[j] + beta2[1] * (1 - grp[j]), 
                            (1 - rho_eps^2) * tau_eps)
      
      for(i in 2:N_years){
        epsilon1[i,j] ~ dnorm(delta1[j] + beta1[i] * grp[j] + beta2[i] * (1 - grp[j]) + 
                              rho_eps * epsilon1[i-1, j], tau_eps)
      }
    }
    
    # Priors for group selection
    for(j in 1:N_sts){
      logit_pi[j] ~ dnorm(0, 1e-4)  # Prior on logit(pi)
      pi[j] <- exp(logit_pi[j]) / (exp(logit_pi[j]) + 1)  # Inverse logit
      grp[j] ~ dbern(pi[j])  # Group assignment - 0 or 1
    }
        
    # Priors
    alpha1 ~ dnorm(0, 1e-4)
    alpha2 ~ dnorm(0, 1e-4)
    tau_global ~ dgamma(0.01, 0.01)

    rho_beta1 ~ dunif(-1, 1)   # Uniform prior for rho_beta1 -- global AR(1) for beta1
    rho_beta2 ~ dunif(-1, 1)   # Uniform prior for rho_beta2 -- global AR(1) for beta2
    rho_eps ~ dunif(-1, 1)     # Prior for rho_eps same for all STs
    
    tau_beta1 ~ dgamma(3, 2)   # Tight prior for tau for beta1
    tau_beta2 ~ dgamma(3, 2)   # Tight prior for tau for beta2
    tau_eps ~ dgamma(3, 2)     # Tight prior for tau, shared for all serotypes
    
    for(j in 1:N_sts){
      delta1[j] ~ dnorm(0, tau_global)  # Serotype means centered around 0
      r[j] ~ dunif(0, 250)  # Serotype dispersion parameter
    }
  }
"

Create JAGS model object

mod9 <- jags.model(textConnection(mod9_string),
                   data = list("N_IPD" = mod9_mat,
                               "N_years" = max(mod9_df$year_n),
                               "N_sts" = ncol(mod9_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

params9 <- c("alpha1", "epsilon1", "delta1", "beta1", "beta2", "tau_beta1", "lambda", "grp")

mod9_postsamp <- coda.samples(mod9, params9, n.iter = 50000)

# save(mod9_postsamp, file = "data/interim/mod9_postsamp.rda")
load("data/interim/mod9_postsamp.rda")

beta_vars9 <- grep("^beta", varnames(mod9_postsamp), value = TRUE)
for (i in seq_along(beta_vars9)) {
  cat(paste0("#### ",i," \n"))
  traceplot(mod9_postsamp[, beta_vars9[i]])
  cat("```\n")
}

1

#### 2 <img src="analysis2_files/figure-html/unnamed-chunk-75-2.png" width="672" /> #### 3 #### 4 <img src="analysis2_files/figure-html/unnamed-chunk-75-4.png" width="672" /> #### 5 #### 6 <img src="analysis2_files/figure-html/unnamed-chunk-75-6.png" width="672" /> #### 7 #### 8 <img src="analysis2_files/figure-html/unnamed-chunk-75-8.png" width="672" /> #### 9 #### 10 <img src="analysis2_files/figure-html/unnamed-chunk-75-10.png" width="672" /> #### 11 #### 12 <img src="analysis2_files/figure-html/unnamed-chunk-75-12.png" width="672" /> #### 13 #### 14 <img src="analysis2_files/figure-html/unnamed-chunk-75-14.png" width="672" /> #### 15 #### 16 <img src="analysis2_files/figure-html/unnamed-chunk-75-16.png" width="672" /> #### 17 #### 18 <img src="analysis2_files/figure-html/unnamed-chunk-75-18.png" width="672" /> #### 19 #### 20 <img src="analysis2_files/figure-html/unnamed-chunk-75-20.png" width="672" /> #### 21 #### 22 <img src="analysis2_files/figure-html/unnamed-chunk-75-22.png" width="672" /> #### 23 #### 24 <img src="analysis2_files/figure-html/unnamed-chunk-75-24.png" width="672" /> #### 25 #### 26 <img src="analysis2_files/figure-html/unnamed-chunk-75-26.png" width="672" /> #### 27 #### 28 <img src="analysis2_files/figure-html/unnamed-chunk-75-28.png" width="672" /> #### 29 #### 30 <img src="analysis2_files/figure-html/unnamed-chunk-75-30.png" width="672" /> #### 31 #### 32 <img src="analysis2_files/figure-html/unnamed-chunk-75-32.png" width="672" /> #### 33 #### 34 <img src="analysis2_files/figure-html/unnamed-chunk-75-34.png" width="672" /> #### 35 #### 36 <img src="analysis2_files/figure-html/unnamed-chunk-75-36.png" width="672" /> #### 37 #### 38 <img src="analysis2_files/figure-html/unnamed-chunk-75-38.png" width="672" /> #### 39 #### 40 <img src="analysis2_files/figure-html/unnamed-chunk-75-40.png" width="672" /> #### 41 #### 42 <img src="analysis2_files/figure-html/unnamed-chunk-75-42.png" width="672" /> #### 43 #### 44 <img src="analysis2_files/figure-html/unnamed-chunk-75-44.png" width="672" /> #### 45 #### 46 <img src="analysis2_files/figure-html/unnamed-chunk-75-46.png" width="672" /> #### 47 #### 48 <img src="analysis2_files/figure-html/unnamed-chunk-75-48.png" width="672" /> #### 49 #### 50 <img src="analysis2_files/figure-html/unnamed-chunk-75-50.png" width="672" />

AR1 groups

st_mapping9 <- setNames(colnames(mod9_mat), 1:length(pcv13_st))
year_mapping9 <- setNames(unique(mod9_df$year), 1:length(unique(mod9_df$year)))

groups9 <- gather_draws(mod9_postsamp, grp[j]) |> 
  rename(st = j) |> 
  mutate(st = st_mapping9[as.character(st)])

groups9_summary <- groups9 |> 
  summarise(mean_grp = mean(.value))  |> 
  mutate(assigned_group = ifelse(mean_grp > 0.5, "Group 1 (beta1)", "Group 2 (beta2)"))
## `summarise()` has grouped output by 'st'. You can override using the `.groups`
## argument.
groups9_summary |> 
  group_by(assigned_group) |> 
  summarise(serotypes = paste0(st, collapse = ", "))

Plot obs vs pred

mod9_summary <- gather_draws(mod9_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping9[as.character(st)],
         year = year_mapping9[as.character(year)]) |> 
  left_join(mod9_df, by = c("year", "st"))

A

mod9_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

B

mod9_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

C

mod9_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod9_postsamp, ask = TRUE)

Plot trajectory of global AR(1)s

AR1 - beta1

mod9_beta1 <- gather_draws(mod9_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping9[as.character(i)])

mod9_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1) - Beta 1") +
  theme_minimal()

AR1 - beta2

mod9_beta2 <- gather_draws(mod9_postsamp, beta2[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping9[as.character(i)])

mod9_beta2 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1) - Beta 2") +
  theme_minimal()

10: Mixture model with PCV-7

(did not converge)

Create JAGS model object

mod10 <- jags.model(textConnection(mod9_string),
                   data = list("N_IPD" = mod1_mat,
                               "N_years" = max(mod1_df$year_n),
                               "N_sts" = ncol(mod1_mat)),
                   inits = list(inits1, inits2, inits3),
                   n.adapt = 10000,
                   n.chains = 3)

Sample posteriors

mod10_postsamp <- coda.samples(mod10, params9, n.iter = 50000)

# save(mod10_postsamp, file = "data/interim/mod10_postsamp.rda")
load("data/interim/mod10_postsamp.rda")

beta_vars10 <- grep("^beta", varnames(mod10_postsamp), value = TRUE)
for (i in seq_along(beta_vars10)) {
  cat(paste0("#### ",i," \n"))
  traceplot(mod10_postsamp[, beta_vars10[i]])
  cat("```\n")
}

1

#### 2 <img src="analysis2_files/figure-html/unnamed-chunk-87-2.png" width="672" /> #### 3 #### 4 <img src="analysis2_files/figure-html/unnamed-chunk-87-4.png" width="672" /> #### 5 #### 6 <img src="analysis2_files/figure-html/unnamed-chunk-87-6.png" width="672" /> #### 7 #### 8 <img src="analysis2_files/figure-html/unnamed-chunk-87-8.png" width="672" /> #### 9 #### 10 <img src="analysis2_files/figure-html/unnamed-chunk-87-10.png" width="672" /> #### 11 #### 12 <img src="analysis2_files/figure-html/unnamed-chunk-87-12.png" width="672" /> #### 13 #### 14 <img src="analysis2_files/figure-html/unnamed-chunk-87-14.png" width="672" /> #### 15 #### 16 <img src="analysis2_files/figure-html/unnamed-chunk-87-16.png" width="672" /> #### 17 #### 18 <img src="analysis2_files/figure-html/unnamed-chunk-87-18.png" width="672" /> #### 19 #### 20 <img src="analysis2_files/figure-html/unnamed-chunk-87-20.png" width="672" /> #### 21 #### 22 <img src="analysis2_files/figure-html/unnamed-chunk-87-22.png" width="672" /> #### 23 #### 24 <img src="analysis2_files/figure-html/unnamed-chunk-87-24.png" width="672" /> #### 25 #### 26 <img src="analysis2_files/figure-html/unnamed-chunk-87-26.png" width="672" /> #### 27 #### 28 <img src="analysis2_files/figure-html/unnamed-chunk-87-28.png" width="672" /> #### 29 #### 30 <img src="analysis2_files/figure-html/unnamed-chunk-87-30.png" width="672" /> #### 31 #### 32 <img src="analysis2_files/figure-html/unnamed-chunk-87-32.png" width="672" /> #### 33 #### 34 <img src="analysis2_files/figure-html/unnamed-chunk-87-34.png" width="672" /> #### 35 #### 36 <img src="analysis2_files/figure-html/unnamed-chunk-87-36.png" width="672" /> #### 37 #### 38 <img src="analysis2_files/figure-html/unnamed-chunk-87-38.png" width="672" /> #### 39 #### 40 <img src="analysis2_files/figure-html/unnamed-chunk-87-40.png" width="672" /> #### 41 #### 42 <img src="analysis2_files/figure-html/unnamed-chunk-87-42.png" width="672" /> #### 43 #### 44 <img src="analysis2_files/figure-html/unnamed-chunk-87-44.png" width="672" /> #### 45 #### 46 <img src="analysis2_files/figure-html/unnamed-chunk-87-46.png" width="672" /> #### 47 #### 48 <img src="analysis2_files/figure-html/unnamed-chunk-87-48.png" width="672" /> #### 49 #### 50 <img src="analysis2_files/figure-html/unnamed-chunk-87-50.png" width="672" />

AR1 groups

groups10 <- gather_draws(mod10_postsamp, grp[j]) |> 
  rename(st = j) |> 
  mutate(st = st_mapping[as.character(st)])

groups10_summary <- groups10 |> 
  summarise(mean_grp = mean(.value))  |> 
  mutate(assigned_group = ifelse(mean_grp > 0.5, "Group 1 (beta1)", "Group 2 (beta2)"))
## `summarise()` has grouped output by 'st'. You can override using the `.groups`
## argument.
groups10_summary |> 
  group_by(assigned_group) |> 
  summarise(serotypes = paste0(st, collapse = ", "))

Plot obs vs pred

mod10_summary <- gather_draws(mod10_postsamp, lambda[i,j]) |> 
  median_hdi() |> 
  rename(year = i,
         st = j) |> 
  mutate(st = st_mapping[as.character(st)],
         year = year_mapping[as.character(year)]) |> 
  left_join(mod1_df, by = c("year", "st"))

A

mod10_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

B

mod10_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st) +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal() +
  coord_trans(y = "log1p")

C

mod10_summary |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_point(aes(y = N_IPD, color = "Observed")) +
  geom_line(aes(color = "Fitted")) +
  facet_wrap(~st, scales = "free_y") +
  labs(x = "", y = "# IPD", color = "") +
  theme_minimal()

# run in console
# plot(mod9_postsamp, ask = TRUE)

Plot trajectory of global AR(1)s

AR1 - beta1

mod10_beta1 <- gather_draws(mod10_postsamp, beta1[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping[as.character(i)])

mod10_beta1 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1) - Beta 1") +
  theme_minimal()

AR1 - beta2

mod10_beta2 <- gather_draws(mod10_postsamp, beta2[i]) |> 
  median_hdi() |> 
  mutate(year = year_mapping[as.character(i)])

mod10_beta2 |> 
  ggplot(aes(x = year, y = .value)) +
  geom_ribbon(aes(x = year, ymin = .lower, ymax = .upper), alpha = 0.2) +
  geom_line() +
  labs(x = "", y = "Global AR(1) - Beta 2") +
  theme_minimal()